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ABSTRACT 

Current plans call for the first Terrestrial Planet Finder mission, TPF-C, to 
be a monolithic space telescope with a coronagraph for achieving high contrast. 
The coronagraph removes the diffracted starlight allowing the nearby planet to be 
detected. In this paper, we present a model of the planet measurement and noise 
statistics. We utilize this model to develop two planet detection algorithms, one 
based on matched filtering of the PSF and one using Bayesian techniques. These 
models are used to formulate integration time estimates for a planet detection 
with desired small probabilities of false alarms and missed detections. 

Subject headings: planet finding, bayesian, hypothesis testing, completeness, 
TPF 



1. Introduction 

In 2004, NASA announced that it plans to launch the Terrestrial Planet Finder-C (TPF- 
C), a space based visible light telescope equipped with a coronagraph to detect earthlike 
planets around other stars, by 2016. Work is proceeding on developing a baseline design 
that would meet the scientific requirements. This currently consists of a 8 x 3.5 meter 
off-axis Cassegrain optical telescope combined with a starlight suppression system (SSS) to 
achieve the needed starlight rejection. The SSS is equipped to utilize either conventional Lyot 
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coronagraphs or pupil masks (shaped or apodized). At Princeton, we have been studying 
what we call shaped pupil coronagraphs as an alternative to the traditional Lyot design. 
(Kasdin et al. 2003; Vanderbei et al. 2003a,b; Kasdin et al. 2005; Vanderbei et al. 2004) 
These are coronagraphs with apodized entrance pupils that rely solely on one/zero binary 
openings, resulting in a more manufacturable and robust design. 

Critical to TPF-C's success is the ability to survey as many stars as possible and max- 
imize the chances of detecting potential planets. To do this, a thorough understanding of 
the data analysis is required in order to obtain reasonable estimates of the integration time 
for each observation. Brown (2005) discusses at length the question of survey completeness 
and its relation to integration time. In this paper, we focus on a more careful and complete 
model of the measurement and noise statistics and develop a probabilistic detection criteria. 
We also analyze a Bayesian detection approach to see if detection time improvements can 
be made. These expressions can be used for further completeness studies and estimates of 
the efficacy of the TPF-C search program. 



Before describing the detection approach, in this section we develop the needed defi- 
nitions, units, and notation for photometry. We will assume a perfect imaging system and 
scalar Fraunhoffer theory such that the electric field at the image plane of the camera is 
the Fourier Transform of the field entering the camera at the exit pupil of the starlight 
suppression system (SSS), 



where E Q is the electric field at the exit pupil aperture, / is the focal length of the camera, 
A is the wavelength of the light, and the set S represents the exit pupil aperture, which, for 
an open square aperture, as an example, is given by, 



where D is the width of the aperture. 

The optical system prior to the exit pupil of the SSS consists of various reflections and 
correction elements all of which have a certain efficiency If we take Ei(X) to be the uniform 
electric field at the entrance pupil, then E Q can be written, 



2. Point Spread Function 




(1) 



S = {(x,y): -D/2 <x< D/2, -D/2 <y< D/2} 



E (X,Y,\) = r ] A(X,Y)E i (X) 
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where 77 < 1 is the overall attenuation of the optical elements up to the exit pupil (including 
any reduction due to an image plane mask), A(X,Y) is the apodization of the exit pupil, 
and we have explicitly shown the dependence on wavelength. For a classical or bandlimited 
coronagraph, A(X, Y) is zero-one valued and just the Lyot stop shape. For a shaped pupil 
coronagraph, A(X, Y) represents the entire starlight suppression and is also zero-one valued. 

For a point source, the image plane detector measures the irradiance of the arriving field, 
If = T^\Ef\ 2 , in photon sec -1 m -2 [imr 1 where c denotes the speed of light, h is Planck's 
constant and e Q is the permittivity of free space. Squaring the magnitude of the electric field 
in Eq. (1) gives the irradiance at the detector, 

I f (X,u,v) = V ^\E t (\)fP(\,u,v) 

= rfl^P^v) (2) 

where = ^\Ei{\)\ 2 is the irradiance of the point source at the reference wavelength 

and P(u,v) is the dimensionless point spread function (PSF) associated with the exit pupil 
of the SSS, 

2 

(3) 

We note that the value of the integral in Eq. (3) at the origin is the square of the integral 
of the apodization, which for an open aperture (or shaped pupil or Lyot stop) is the square 
of the area of the exit pupil. We use A to represent the area of the entrance pupil of the 
telescope and A c to represent the area of the exit pupil of the coronagraph (Lyot stop or 
shaped pupil). For an open telescope with no coronagraph, A c = A. For apodized pupils, 
we call A c the pseudo-area as it includes the reduction in throughput due to the continuous 
apodization. 

There are three expected sources that produce electric field entering the telescope, the 
parent star (I s = ^\E S (X)\ 2 ), the planet (I p = t^|-E p (A)| 2 ) and the exo- zodiacal dust. While 
the first two sources are point sources, the exo-zodi is a uniform extended source with flux 
jF e2 (A) in units of photon sec -1 m~ 2 /an" 1 ster -1 . 

For a planet offset on the sky by the angles (0 x ,6 y ) the electric field at the entrance 
pupil becomes E p (\)e 27T ^ 9xX ^ x+eyy ^ x ^ and by the Fourier Shift Theorem the irradiance at the 
detector becomes, 

I f (\, u, v) = r, 2 I p {\)P{\ u - 9J, v - 9 y f) (4) 

For a diffraction limited system, we can assume that the uniform exozodical flux also 
creates a uniform flux density on the image plane and the details of the PSF are unimportant 
(except for perhaps a slight broadening at the edges). Thus, the irradiance at a single pixel, 



P(X,u,v) 



\ 2 P 



J ^A(x,y)e-^ {xu+yv) dxdy 



-4- 



is given by I eZij = jF e2 (A)f2ijAo; i:; -, where Qij is the solid angle on the sky subtended 
by a single pixel and Aa^ is the area of a pixel. 



In this section we develop a model for the detector measurements and corresponding 
noise statistics. This is essential for developing the detection algorithm and corresponding 
detection probabilities. We assume here that the wavefront control system has successfully 
reduced the aberrated light to a uniform background in the discovery region or that we 
are dominated by the exozodiacal background. Thus, the intensity at the planet location 
consists of the planet PSF from Eq. 4 plus a uniform background, If, = I ss + I sc + I ez + I r , 
where I ss is due to the residual suppressed starlight, I sc is due to scattered light, I ez is the 
exo-zodiacal dust background, and I r is the effective readout noise. We assume that because 
of cosmic ray concerns the eventual signal will consist of a stack of many frames read out at 
a faster rate and thus it is sensible to treat the readout noise as absorbed into the uniform 
background. We also assume that the readout noise is small compared to the other sources 
of background noise. In fact, it is expected that the dominant background source will be the 
exo-zodiacal light, perhaps as much as 3 to 10 times the brightness of the planet. 

The photon arrival at each pixel is determined by a Poisson distribution. That is, the 
probability of a given measurement, at pixel and time t is given by: 



where is the mean arrival rate of the photons, or intensity, at pixel in photons/sec. 
Thus, if z(t) is the actual measured number of counts at any pixel, 



where I is the intensity at that pixel. 

The photon arrival rate, or intensity, in photons/sec, at pixel is the integrated 

irradiance over the pixel area and a reference passband about a reference wavelength A r , 
which for the planet is given by, 



3. 



Measurements 



Pr[Z = 




£{z(t)} 
£{{z{t) - pt(t)) 2 } 



It = fx(t) 

It = a 2 (t) 
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where Aa^ represents the region covered by pixel and AA is the reference wavelength 
band given by the largest band of interest. 1 In general there is not much more that can be 
done without a specific point spread function (and then these integrals can be difficult in 
closed form). However, since we are seeking a simple and conservative metric, we can obtain 
a useful result by assuming that the bandpass is small relative to X r f /D. In this case, the 
pixel intensity reduces to the simpler form, 



la = rj AXI p (X r )Aa 



1 



(AA r /2A r ) 2 



Pa (6) 



where the unitless quantity, Py = f / A P(X r ,u — 9 x f,v — 9 y f)dudv, is the normalized 
integral of the PSF over pixel (i, j), and Aa is the pixel area. Since we are considering small 
wavelength bands, we can safely drop the term in brackets in Eq. (6). 

The final electron count rate at each pixel (also Poisson distributed), including back- 
ground, is then, 

4 = eri 2 A\I p (X r ) AaP(0, 0)Py + erfAXI b Aa (7) 

where e is the quantum efficiency of the camera and we have normalized the PSF by its 
value at (0,0) (so that = Py/P(0, 0) and likewise P(\,u,v) = P(A, u, v)/P(X, 0, 0)). We 
introduce the carat to differentiate between the photon and electron counts. 

Finally, it will become useful later to utilize the contrast ratio, Q, first introduced by 
Brown and Burrows (1990). The contrast ratio is defined as I p P(0 : Q)/h- In words, Q is the 
ratio of the peak of the planetary image (in photon sec^ 1 mr 2 /im" 1 ) to the local surface 
brightness of the background (in the same units). 

It is often helpful, and very common, to non-dimensionalize the image plane coordinates. 
Recalling that P(A r ,0,0) = A 2 c /X 2 .f 2 , the normalized PSF is given by, 



P(u, v) 



P(0, 0) A\ 



J J A{x,y)e- j ^ {xu+yv) dxdy 



(8) 



Next, introducing a reference length, D (the side of a square aperture or the diamter 
of a circular aperture, for instance), defining dimensionless pupil plane coordinates, X = 



1 For TPF-C this widest band of interest, and the one used for this ideal reference case, is planned to be 
from 0.5 to 0.8 micron. Strictly speaking, the variation of the field intensity, I P (X), should be included in the 
integral over wavelength as the Planck curve has significant variations over this wide a band. For simplicity, 
however, we are approximating the response curve as flat here for the purposes of developing the metric. In 
the less idealized case where we measure over narrow bands this approximation becomes much better. 
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x/D,Y — y/D, and a dimensionless angular position in the image plane in units of X/D, 

Du _ Dv 
e " A/' C " Xf 



gives the dimensionless form of the PSF, 



D 4 



c 



J J_A(X,Y)e~^ x ^dXdY 



(9) 



This equation can be modified slightly by introducing the definition sD 2 = A. That is, the 
entrance aperture area can be related to the square of the reference length by the constant 
s. Thus, for example, for square apertures s — 1, for circular apertures s = ir/4, and for 
elliptical apertures s = irr/4, where r is the aspect ratio of the ellipse. The normalized PSF 
in this notation is given by, 



A 2 1 

^•0 = ^ 



c 



J J_A(X,Y)e-^ x ^dXdY 



(10) 



Finally, (w, v) is also replaced in the expression for P i3 - to get the pixel integrated PSF 
in dimensionless units, 

where Aa is the dimensionless pixel area (AcS = ^pAa). Thus, as expected, Pij can be 
computed in either dimensional or non-dimensional coordinates. 

Substituting for P(0, 0) and Aa, the measurement equation can be re-written, 

A 2 

I^t) = e^AXIpiX^Aa^-sPij + e7] 2 AXI b Aa (12) 

Finally, we define the total throughput, T = A c /A, the ratio of the exit pupil pseudo-area 
to entrance pupil area. Substituting into Eq. 12 leaves, 

4(t) = erf AXI p (X r )AaT 2 AsPij + er] 2 AXI b Aa (13) 



Thus, the intensity has been rewritten in terms of the total throughput and dimension- 
less pixel area. 
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4. Planet Detection 

In this section we examine planet detection via a linear combination of pixel measure- 
ments. The usual approach is to do PSF fitting where the irradiance of the planet is estimated 
by fitting a PSF to CCD image using standard photometric approaches (Brown 2005). In 
other words, a planet PSF fit is attempted at every pixel and the results compared to find a 
statistically significant outlier. A linear, unbiased procedure for PSF fitting is described in 
Burrows (2003) and summarized in the appendix. Here, we provide a probabilistic basis for 
using PSF fitting for planet detection and develop unambiguous metrics. We also show that 
PSF fitting is, in fact, the optimal linear detection strategy. 

4.1. PSF Fitting 

As in Burrows (2003), we represent the Poisson counting process by a linear measure- 
ment model, 

z ij = C p P ij + C b + u (14) 

where C p = erfA\I p (\ r )AaT 2 sAt ) C b = er] 2 A\I b Aat, v represents the photon and readout 
noise, and t is the integration time. Thus, C p is approximately the mean photon count at 
the pixel centered on the PSF (since Poo ~ 1) and C b is the mean background photon count 
at any single pixel. 

Given this model in Eq. 14, a variety of methods are available to estimate C p , the 
planet irradiance. The most common approach to photometry is via PSF fitting (also called 
"Matched Filtering"). The standard, linear method is described in Burrows (2003). We 
summarize this derivation in the Appendix and compare to maximum likelihood approaches. 

While Burrows (2003) primarily considers the background limited case, here we must 
assume that the planet signal and the background are of comparable intensity. (It is pos- 
sible that the exo-zodi could be as much as three to four times the intensity of the planet. 
Nevertheless, this is not large enough to justify the background limited regime and thus it 
is best to include both the planet and background noises. Note that Brown (2005) uses the 
background limited results of Burrows (2003) yet includes the planet signal, resulting in an 
inconsistency in his results.) 

While it is possible to consider nonlinear maximum likelihood approaches to PSF fitting, 
these are complicated, provide minimal improvement (if any), and don't have convenient 
expressions for the estimate and variance. We therefore continue to use the linear, unbiased 
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estimator suggested by Burrows (2003), 



a _ T.ij{ z ij -c b )Pij 

°p - p2 

t—tij ij 



where the sum is taken over some region AS of the image plane covered by the PSF. 



2 



For pixels with a planet, this is an unbiased estimate of the planet signal, C p , with 
variance, 

2 _ C P Tnj Pjj C b 

\/—/ij ij) '-—'ij ij 

For pixels without a planet, the expected value of Eq. 15 is zero and the variance is, 

°l = y^W (17) 
Z^ij ij 

Planet detection is performed by calculating the estimate in Eq. 15 at each pixel in the 
dark hole region. For pixels without planets, the estimate will scatter about zero within 
the range defined by the variance in Eq. 17. For pixels with a planet, and a long enough 
integration time, the planet estimate should exceed this region about zero by some amount 
(see Fig. 4 for an illustrative example). We can plot the estimate for all pixels and look for 
where it exceeds the la bounds defined by Eq. 17 by some factor, usually 5 or greater (a 5a 
detection). 



4.2. Detection Metric and Missed Detection Probability 

Burrows (2003), and later Brown (2005), define the SNR as the ratio of the actual 
planet irradiance, C p (the expected value of the estimate), and the standard deviation of 
the estimate error. While this is a reasonable indicator of the quality of the estimate, it 
says little about how to make a detection decision. The irradiance, C p , is not the signal 



2 Note that throughout we have assumed that the background is uniform in the summation region about 
the PSF (roughly 10 to 20 pixels). We also assume the mean value is known in order to do the subtraction 
in Eq. 15. The uniformity assumption is easily relaxed without changing the results by simply indexing 
Cb, though on this scale the cxo-zodi should be very uniform. However, if the background is unknown the 
problem becomes more complicated. For telescope fixed background (speckle due to errors in the optics), 
the subtraction is accomplished by rotating the telescope and subtracting two measurements, at the expense 
of doubling the background noise. If the zodi is uknown, then a simultaneous estimation is necessary; this 
complicates the results and will be treated in a future paper. 
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nor is the estimate variance the noise. This approach to computing integration time also 
ignores that fact that the estimate of C p is also a random variable and thus has a specific, 
quantifiable probability of failing any detection test. Finally, Brown (2005) uses both the 
planet and background photon noise in the definition of the SNR. As just described, only 
the background is needed to perform a detection. 

For the detection problem, we are interested in a metric that can be used to decide 
whether a planet exists in a given region. This can be used to determine a maximum 
integration time before moving on. The signal in this case is the estimate of the irradiance, 
C p (from Eq. 15), and the noise is given by the background standard deviation on pixels 
without a planet (from Eq. 17). Thus, a detection occurs when the estimate exceeds the 
background by some amount. The resulting signal-to-noise ratio metric, C p /<7 b , is a random 
variable and we can therefore define a confidence interval for the detection decision. We 
choose an integration time such that the probability that the SNR, assuming a planet 
exists, exceeds a given threshold, K, is equal to a desired value, P, 



Pr |S > ^|a planet j = P (18) 

For example, for 99% confidence, P = 0.99. That is, if the correct integration time is used 
and no planet is seen, then we can say with 99% confidence that there is no planet at that 
location (or at least, that there is no planet within Ka b of the background noise). 

To compute this probability, we first use Eqs. 15 and 17, to find the SNR, 

C p _ J2ij( z ij - Cb)Pi 



(19) 



Unfortunately, the probability mass function of C v j a b is extremely difficult to compute 
using a Poisson process for the Zij. However, for the integration times being considered, the 
arrival rate is high enough that the Poisson process is well approximated by a Normal one and 
this assumption allows for a straightforward, and adequate, estimate of the integration time. 
Also, by the Central Limit Theorem, the sum in Eq. 19 normalizes the random variable. 
From Eq. 19 the mean and variance are easily computed, 



PLSNR = \/C p Qj2 P Z ( 20 ) 

a SNR = 1 + sp p2 ( 21 ) 
2-1 ij 

where we have used the definition of Q introduced earlier. 
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This allows the probability in Eq. 18 to be rewritten in terms of a normalized variable, 

Pr [z > K ~ fISNR \ & planet) = P (22) 

where Z is a standard normal random variable with zero mean and unit variance. The 
detection criteria is then given by the probability integral, 

1-<S>( K -» SNR ) = P (23) 

V ° SNR J 

or 

$ |^-_/W) = 1 _p = p MD (24) 



&SNR 



where Pmd = 1 — P is the probability of a missed detection and &(z) is the gaussian 
distribution function, 

= Pr[Z < z] = J Z ^e^dy = I (l + ed(z/y/2)) (25) 



In other words, the integration time is chosen so that — flSNR reaches a value, 7, such 

' & &SNR ' ' ' 

that $(7) = 1 — P. For example, a 0.1% missed detection rate (Pmd = 0.001) corresponds 
to 7 = —3.1. The detection criteria can now be written, 

= 7 (26) 



1 + 



To solve for the needed integration time to achieve the desired false alarm and missed detec- 
tion rates, this equation is first solved for C p (since Q is independent of t) and then divided 
through to find the integration time, 



t = : ^ 1 = 7 (07) 

erj 2 A\I p (\ r )T 2 As QAa^Pl 



For consistency with Brown and Burrows (1990) and Burrows (2003), and for computa- 
tional convenience, we make one more modification by introducing the sharpness. Sharpness 
is defined in Burrows (2003) as ^ = . P? / '(^ Pij) 2 . Using sharpness has the value of 
making the calculation less sensitive to the location of the PSF within a pixel. Sharpness 



- 11 - 



is a unique characteristic of the PSF of any coronagraph. With this substitution, Eq. 27 
becomes, 



1 ( K - 7\/l + f ) 



* = 7;- — (28) 

v p. 3 . 

where we have also defined (3 = er] 2 AXI p (X r )TA, Q = Q^^Pij, 2 = ^ ^ 3 and = 

TsAa^Pjj = Ts J f AS P(£X)d£,d(. Q is the ratio of the total photon count from the 
planet over the integration area, AS*, to the number of background counts in a single pixel. 
Ta is what we called the Airy Throughput in Vanderbei et al. (2003b). It is the throughput 
times the integral of the normalized PSF over the integration area AS*. It is straightforward 
to show, using Parseval's Theorem, that s f J As P '(£, C) d£d( is one f° r an y open aperture 
assuming that AS covers the entire image plane. For the expected smaller integration area, 
Ta is always less than T. It is thus the fraction of the total throughput achieved for the 
given region of the PSF used in detection. This difference can be significant. For an open, 
circular aperture, the Airy throughput corresponding to the first null is only 84.5%. If only 
the FWHM is used, then this drops to 53%, almost a factor of 2 in integration time. For 
shaped pupils, this is also important since most shaped pupil PSFs move a significant portion 
of the light outside the outer working angle; this is light not used in the planet detection. 

This trend of the Airy Throughput would imply that it is always best to use as many 
pixels as possible to cover as much of the planet's PSF as possible. However, the remaining 
terms that depend upon Q, and S have the opposite trend, suggesting that a smaller 
number of pixels around the central core is better (since more pixels tend to include more 
background light but less planet light because of the drop off of the PSF). The result is 
that beyond a certain region of the core there is no improvement in integration time. This 
region size depends upon the specific PSF and the size of the pixels. For an Airy PSF it is 
roughly 0.75 X/D. Figure 1 shows two plots of the normalized integration time ((3tT) versus 
width of the core used in the detection for four different pixel sizes and a background limited 
detection, where Q = 1/3. While there is an advantage to oversampling the PSF, the gain 
is slight. For smaller background, and larger Q, slightly less of the core is necessary but the 
difference in integration time is also small. 

This also raises the question of the location of the PSF within the pixel. Up to now, this 
has been ignored. Figure 1 shows that by using sharpness the integration time expression is 
reasonably insensitive to the PSF location relative to the pixel center. However, in practice, 
the estimate of C v in Eq. 15 assumes that Pij is known at every pixel, introducing an 
ambiguity. One solution is to attempt a fit for many values of P^ as the PSF moves relative 
to a pixel center. Alternatively, we could compute the sensitivity of the detection to errors in 
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Pij and increase the integration time accordingly. Finally, we could treat the actual value of 
Pij as a random variable with uniform distribution between its value for the PSF landing on 
a pixel center and its value for the PSF landing on a pixel edge. This probability distribution 
can then be included in the probability calculation in Eq. 18. We defer this approach to 
future work. 



4.3. Choice of K and False Alarm Probability 

Equation 28 provides the integration time to achieve a desired small missed detection 
rate for a given detection threshold, K. However, the choice of K still seems arbitrary. Most 
astronomical applications select K between 5 and 10 (while ignoring the missed detection 
probability embedded in 7). In reality, K has a probabilistic basis as well; it is a measure of 
the probability of false alarms. In fact, the best approach is to pick the integration time, t, 
and the threshold, K, to achieve a desired balance between the two types of errors. 

The probability of false alarms is given by an expression similar to that in Eq. 18, 

Pr |^ > K I no planet J = P FA (29) 

The expression for C p /<Jb is the same as that in Eq. 19, but since there is no planet the mean 
and variance differ, 

UFA = (30) 
a\ A = 1 (31) 

The SNR when there is no planet in a pixel is, in fact, well approximated by a standard 
normal random variable. The false alarm probability is then given simply by, 

$(X) = 1 - P FA (32) 

Thus, Eq. 32 is used to select the threshold, K, corresponding to a desired probability of false 
alarms and Eq. 28 is then used to select the integration time to achieve a desired probability 
of missed detections. 

Since we expect at best to only find one planet per system surveyed, most pixels will not 
have a planet. This places a high premium on keeping the false alarm rate low (a large K), 
otherwise much time would be wasted on follow up of false planets. Nevertheless, a K greater 
than 5 results in a false alarm probability lower than necessary (e.g., <£>(5) = 0.9999997, a 
false alarm probability of 3e-7). The result is a potentially longer integration time than 



13 



necessary, dramatically affecting the mission. A reasonable approach is to select the false 
alarm probability such that there is less than one false alarm per observed system. For 
instance, for a TPF mission that surveys 100 systems, a total false alarm rate of 1% would 
result in 1 false alarm over the course of the mission. Since each planet detection will be 
followed by deep surveys for characterization, the false alarm would be quickly discovered 
and have minimal impact on mission life. If we assume on the order of 3000 measurements 
in a given system, we therefore require that P F a x 3000 ~ 0.1 or P FA ~ 3 x 10~ 5 . This 
corresponds to K — 4. 

The missed detection probability is critical for a mission with so few surveyed systems 
and unknown planet occurrence rate, 77®. It should be made low enough to avoid a single 
missed detection over the course of the mission. A missed detection probability of less than 
0.1% is most likely needed, requiring 7 = —3.1. 

With these numbers, the integration time equation becomes, 



There are a number of notable differences between this expression for integration time 
and that in Brown (2005). First, Brown (2005) selected the SNR threshold by implicitly only 
considering the false alarm probability. He ignored the missed detection probability (7) and 
the potentially significant contribution to the integration time in Eq. 28. (Physically, this 
corresponds to basing the detection decision only on the mean value of the SNR which leaves 
a significant probability that a particular SNR may fall below the threshold. Mathematically, 
this means 7 = which corresponds to a missed detection probability of 50%.) Brown (2005) 
also neglected to account for the possible reduction in Airy throughput due to only using 
the center of the PSF. 

Another difference is that Brown (2005) included the planet photon noise in the SNR 
expression. We describe above why only the background noise should be used in the SNR 
definition. 

Also note that the needed integration time depends upon the actual planet irradiance 
(in f3). Strictly speaking, this is itself a random variable as the planet irradiance could vary 
quite widely among systems. While it is possible to include a density function for I p and 
marginalize the SNR, this makes the problem substantially more complicated. It is also 
not clear that an a priori density for the planet irradiance can easily be found. A more 
reasonable approach is to use a worst case irradiance in Eq. 28 for the most conservative 
integration time. This will become important again in the Bayesian analysis later. 




(33) 
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4.4. Numerical Example 

In this section we present a simple numerical calculation of the integration time for a 
typical system. For simplicity, we ignore the ellipticity of the TPF-C entrance aperture. To 
begin, we also ignore the coronagraph (and any affect on the PSF) and consider the shortest 
possible integration time to detect a planet given the telescope aperture. In other words, 
A c = A , the entrance aperture area. I assume the entrance area, A, is the same as the 8m 
x 3.5m TPF-C design (22 m 2 ). Assuming a 5th magnitude star and a planet with Amag 
of 25, the irradiance is approximately I p = 9.5 x 10~ 9 photons cm -2 nm" 1 sec" 1 . With a 
QE of 0.8 AA = 100 nm, and if = 0.33, (3 = 0.06 photons/sec. The integration time then 
becomes, 

(4 + 3.1^1 + f 

t 18^ = 

QT A * 

Assuming a fairly conservative value of Q = 1/3 (large exo-zodi), we can use Fig. 1 to 
obtain the normalized integration of roughly 400 for a critically sampled Airy function. The 
resulting integration time becomes 7200 seconds or 2 hours. 

A more realistic case would include the effect of a Lyot stop (or shaped pupil). The 
Lyot stop reduces T and Ta (by 0.3 for the current design). If we assume that a detector 
with larger pixels is used so that the Airy function is still critically sampled, the sharpness, 
S, and Q are all unchanged. The result is an increase in the integration time by a factor 
of 1/0. 3 2 or about 11. This implies a required integration time of 22 hours. 



4.5. General Linear Detection & Optimality 

The specific irradiance estimate in Eq. 15 can be easily generalized into a form that 
covers all possible linear processing techniques, 

Z = J^W lJ (z tJ -C b ) (35) 

ij 

for some unknown set of weights, Wij, over a region AS* of the image. Following the same 
process as before, the signal-to-noise ratio is, 

— = — / (36) 
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with mean /i S nr = y/C p Q W^P^/ W ij and variance a 2 SNR = 1 + ^V-'^''- % 
the same probabilistic argument, the detection condition becomes, 



i + 



7 (37) 



As before, this can be solved for C P and then t, 



(3 QAaiZWijPij) 2 



(38) 



One observation from Eq. 38 is that the weights can be multiplied by an arbitrary 
constant without effecting the integration time (and thus the detection). This is because 
any constant weight divides out when forming the SNR statistic. Thus, while the constant 
term Pf- in Eq. 15 is important for an unbiased estimate of C p (photometry), it is not 
necessary for a simple detection. 

Equation 38 can be used to find the optimal weights resulting in a minimum integration 
time. These turn out to be very close to PSF fitting in the last section. 

This minimum is found by differentiating Eq. 38 with respect to and setting equal to 
zero. Unfortunately, the exact problem is difficult to solve because of the term in the radical. 
However, for Q <C 1, this term is small and the integration time is very well approximated 
by, 

1 (^-7) 2 E^- , , 

PQAa&WijPij) 2 1 ] 

The resulting minimum is found from, 



d ( £wg 



(40) 



Completing the differentiation results in the condition for W^, 

' £ mAj 



Wij = Pn== — ^- (41) 



While difficult to solve in general, it is easy to show that this equation is satisfied 
by the PSF fitting weights used in the previous section, = Pij. Thus, for small Q 
(large background), the PSF fitting approach is the approximately optimal linear detection 
technique. Any improvement from searching for the exact optimal weights is only very slight. 
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5. Bayesian Detection 

We showed above that PSF fitting for detection was the optimal, that is, shortest time, 
detection among all linear algorithms. It can also be shown that under certain assumptions 
it achieves the Cramer-Rao bound and is thus the best possible estimate of the irraadiance. 
However, for the detection problem we are not necessarily interested in planetary photom- 
etry, just a determination of whether a planet exists. PSF Fitting may then provide more 
information than is necessary. It is reasonable, therefore to question whether there might be 
faster algorithms for making the detection decision. In this section we explore an alternative 
hypothesis testing approach using Bayesian methods using the same measurement model as 
in Section 3 in an effort to answer this question. 



5.1. Detection Criteron 

To determine whether a potential planet is present in a set of pixels AS, centered at 
pixel n, we compare the probabilities of the following two alternate hypotheses: 

Hq : there is no planet in the n — th pixel, 
if" : there is a planet in the n — th pixel. 

Each hypothesis, Hi, has the associated model Mf 

Mq : Xij = Ch 

: = CpPij + Ch 

where X^ is the mean arrival rate of photons in the Poisson distribution at pixel located 
in AS. 

We now compare the posterior probabilities of the two models conditioned on the col- 
lected data and decide that there is a planet in pixel n if the odds ratio, O™ , favors H™ over 
H n : 

= p(M?\Z) = p(M?)p(Z\M?) = p(Z\M?) 
10 p(MS\Z) p(MZ)p(Z\MZ) p(Z\MZ) { > 

where Z is the set of collected measurements and we have used Bayes' rule to write the odds 
ratio as the product of the ratio of prior probabilities, R = p(M") /p(Mq), and the ratio 
of likelihoods, p(Z\M™) /p(Z\Mq). The dependence on the ratio of prior probabilities of 
there being a planet is troubling at first since it appears that the effectiveness of the test can 
be greatly affected by this essentially unknown quantity (alternatively, one can view this as 
perhaps an attractive feature as prior information on the frequency becomes available). We 
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will show later that this dependence is easily removed by balancing the sensitivity between 
missed detections and false alarms. 

One possible and simple decision criterion would be that O™ > 1 favors if™, thus 
indicating that there is a planet whose true location (* is inside AS. However, a better 
approach to the decision is to base it on a different value of the ratio derived from particular 
weighting criteria. This is called a "Loss-based decision". The specific value of the Odds 
Ratio for the decision is selected to balance the relative priorities between missed detections 
and false alarms (note the similarity to the PSF fitting test above). In most applications 
this process is used to optimize the power of the test (via the Receiver Operator Curve or 
ROC). Here, however, as discussed above, there is a large incentive to reduce the probability 
of false alarms well below that of missed detections rather than balance the two types of 
errors. In what follows, analytical expressions are developed to quantitatively trade off these 
two probabilities. 



5.2. Likelihoods and Priors 

In general, the probability of a given data set depends on various physical parameters of 
the system. For instance, the mean photon rate in the Poisson distribution depends upon the 
irradiance of the planet and its location in the pixel. As mentioned in the previous section, 
this is, strictly speaking, a random variable. The usual approach to computing the likelihood 
functions is to marginalize the probabilities to remove the parameter dependence altogether. 
However, given our lack of knowledge regarding the planet distribution, and that our goal 
herein is to simply find a conservative estimate of the needed integration time for detection, 
it is more sensible to begin with a worst case analysis. Thus, we assume the parameters 
(such as planet irradiance) are known and take on their most pessimistic value (for instance, 
for planets below a certain limiting Amag we assume a detection is impossible). If an actual 
planet turns out to be brighter, then we have simply integrated longer than necessary, but 
we will not miss the detection. In short, we trade conservatism for our lack of knowledge of 
the specific distribution of the priors. 

The odds ratio now becomes the ratio of the likelihood functions times the prior prob- 
abilities, 

O n w = (43) 
where the likelihood functions are given by L\ = p(Z\M™) and L = p(Z\Mq). Each of these 
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probability distributions per pixel are Poisson with parameter A^-, 

L =U A-V?' X ( 44 ) 

IJ J 

Thus, for a set of data per pixel, Zij, we form the Odds Ratio using our worst case value for 
the planet irradiance, C p , and the contrast ratio, Q, where I have used the bar to contrast 
the value used in the detection algorithm from the true values, C p and Q. The background 
irradiance, Cb, is still assumed to be known (and uniform). Substituting into Eq. 43, the 
Odds ratio becomes, 

O™ = Re-^ AP + OI' :j y- (45) 

ij 

where AP = Ar ft i s often more convenient to use the log of the Odds Ratio as the 
test statistic, 

logO™ = log R - C P AP + % log(l + QPij) (46) 

ij 

The two constant terms don't affect the statistical test, so it is easier to define a new 
test statistic that removes them, 

X B = £*ylog(l + QPy) (47) 

ij 

Unlike in PSF fitting, not only do we need to know P^, we also need to know Q. This 
begs the question of the sensitivity of the method to errors in our assumed knowledge, Q. 
As mentioned above, for now we use a worst case approach, using the smallest value of I p 
expected (recall that in both of these approaches the background is assumed known). We 
explore the validity of this approach later. 

The test statistic, \i depends upon the measured data and assumed, worst case param- 
eters. What remains is to find a threshold value for x to decide between the two hypotheses. 

5.3. Hypothesis Tests 

Again, the Loss based decision is made based upon the value of the test statistic. That 
is, if x n > 7 then we reject Hq and decide that a planet is present in AS. Otherwise, we 
accept Hq and assume no planet. 7 is selected to balance the two types of errors. 
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5.3.1. Missed Detections 



The probability of a missed detection can be written formally as, 



Pr{ X n < l\M?} 



or 



Pr{J2 % log(l + QPij) <l\M?) = Pi 



MD 



(48) 



To compare to PSF fitting, we look for an expression that allows us to pick an integration 
time to achieve a certain desired probability of missed detection, Pmd- Unfortunately, even 
knowing that Zk is Poisson, forming the probability distribution in Eq. 48 is formidable. 
For a precise result, we must resort to Monte Carlo simulations. However, for the expected 
integration times and arrival rates here, the Poisson distribution is close to Gaussian and, 
as in Section 4, this assumption is useful for obtaining integration time estimates. Thus, 
letting = log(l + QPij), the probability in Eq. 48 can be rewritten in terms of a standard 
normal random variable, 

Pr{Z < = P MD (49) 

0~Z 

where [iz = C p Bij (Pij + 1/Q) and erf = C p Bfj (Pij + 1 /Q) and Z is normally distributed 
with mean and standard deviation 1. Thus, the missed detection equation becomes, 



where, as with PSF fitting, the desired missed detection probability, Pmd, is picked a priori. 
Letting a be the value of the argument such that $(a) = Pmd, the integration time criteria 
becomes, 



In other words, for a desired missed detection probability, the integration time is selected 
such that Eq. 51 is satisfied. 

Substituting for jj,z and oz leaves, 
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7 - C v BijiPij + l/Q) = ayJc^mPii + VQ) 



(52) 



This equation has two unknowns, however, 7 and t. As before, the false alarm rate provides 
a criteria for the threshold, 7, that can be substituted back in to find the integration time. 
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5.3.2. False Alarms 

The probability of a false alarm is similar to Eq. 48, 

Pr{ X n > 1 \M%} = P FA 
or 

Pr{J2 % log(l + QPij) >7|M ™} = Pfa (53) 
As before, this is rewritten in terms of a standard normal random variable, 

Pr{Z' > = 1 - d> ( 1^*) = Pfa (54) 

<jz> V a z< J 

where here \iz' = Yl BijCb and <j\, = Bf^Cb- Letting 5 be the value of the argument such 
that $(5) = 1 - P FA leaves, 

= 5 (55) 

or, 

This equation provides the detection threshold test for x n for a given false alarm probability 
(embodied in the value for 5). 



5.3.3. Integration Time 



The Bayesian integration time can now be found by using the false alarm based expres- 
sion for 7 in Eq. 56 in Eq. 52, 



C„ 



^ - I>>M'j + 1/Q)] = ay/c^BUKj + l/Q) ~ S \fcr^ ( 57 ) 



This can be used to solve for C p as before and subsequently the Bayesian integration time, 



tn = 



eAXI p (X r )TA c s 



QAa (Y.BijPii)' 



(58) 



Finally, we can use the same renormalization as before and write the Bayesian integra- 
tion time in terms of (3, Ta, and Q, 



t B = 



B 



P 



QTa {J2 BijPij/ Bij) ' 



(59) 
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Surprisingly, for background limited systems where Q is small, there is no improvement 
in integration time with the Bayesian approach. This is rather easy to show using Eq. 59. 
For small Q, = log(l + QPij) = QPij- With this subsitution, Eq. 59 becomes exactly 
Eq. 28. Thus, this Bayesian approach, where we assumed we knew the background and used 
a worst case value of I p (i.e., Q), provides no integration time benefit over PSF fitting for 
small Q (large background). Additionally, one would expect a poorer performance (that is, 
a slight increase in the probability of false alarms and missed detections) when Q ^ Q. Also, 
this approach suffers from the same ambiguity regarding the location of the PSF relative to 
a pixel center as PSF fitting, as P^- has been assumed known throughout. 

One of the attractive features of the Bayesian approach is that it provides a methodology 
for handling these unknowns. For example, the uncertainty in Q can be removed by assuming 
a distribution for I p and marginalizing the likelihood function. Likewise, the ambiguity with 
regard to pixel size can also be removed via marginalization. Finally, we can study the 
case where the background is not known by again marginalizing over a distribution for the 
background photon count. An unknown background can also be treated in PSF fitting by 
simultaneously estimating the planet irradiance and background, at the cost of much greater 
complexity and much longer integration times. We defer these comparisons to a future paper. 



6. Monte Carlo Simulation 

As a straightforward verification, we ran 50,000 monte carlo runs for each case, with a 
planet and without. The Poisson arrival rate was computed using the integration time values 
computed in Eqs. 33 and 59. An Airy PSF was used for all simulations and Q was taken 
as 1/3. Figure 2 shows histograms of the results of the signal-to-noise test statistic using 
PSF fitting and Fig. 3 shows the equivalent results for the Bayesian test statistic, x n - Also 
shown on the plots are dotted lines indicating the threshold test levels for the false alarm 
and missed detection probabilities described above. In both cases, there were 40 missed 
detections for the 50,000 runs with planets. For PSF fitting there were 3 false alarms and 
in the Bayesian approach there were 2. These plots confirm the analytical results and also 
justify the Gaussian assumption. 

Figure 4 shows 200 samples each of the PSF fitting test statistic with and without a 
planet. This illustrates how the test statistic indicates the presence of a planet by exceeding 
the zero mean scatter by a threshold amount. 
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Fig. 1. — Normalized integration time versus size of Airy PSF used in detection for various 
pixel sizes. Left PSF centered on single pixel. Right PSF centered on corner of four pixels. 




Fig. 2. — Histograms of 50,000 Poisson simulations of the PSF Fitting signal-to-noise ratio 
test statistic for each case, pixels without a planet (Left) and pixels with a planet (Right). 
Dotted line indicates the threshold test value, K. 
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7. Final Remarks 

In this paper we developed a careful model of the coronagraph measurement and used 
it to define two different detection approaches, matched filtering and Bayesian. For either 
approach, we developed a simple expression for the needed integration time to achieve the 
desired probability of success. These integration time results can now be utilized for mission 
completeness studies. 

The story is far from complete, however. Critical to these studies is an understanding of 
the expected background and the best possible estimates of Q. Additionally, we have ignored 
the question of wavefront error and control. In this simplified model, the background was 
assumed known and uniform. In reality, the background is going to be a mix of exozodiacal 
light and residual wavefront error. The speckle background due to uncorrected wavefront er- 
ror will not be uniform and can potentially confuse the planet detection algorithm. While the 
planet and background can be simultaneously estimated, it complicates the algorithm and 
lengthens the integration time. Also, the additional flexibility of multi-wavelength measure- 
ments has not been explored. Finally, the sensitivity questions of mismatches in assumed 
parameters and PSF alignment on pixels has not yet been studied. These issues are all 
currently being addressed and will be presented in a future paper. Nevertheless, the funda- 
mental model and technique presented here provides an important benchmark to measure 
the potential search capabilities of TPF. 



REFERENCES 

N. J. Kasdin, R. J. Vanderbei, D. N. Spergel, and M. G. Littman. Extrasolar planet find- 
ing via optimal apodized-pupil and shaped-pupil coronagraphs. The Astrophysical 
Journal, 582:1147-1161, January 2003. 

R.J. Vanderbei, D.N. Spergel, and N.J. Kasdin. Spiderweb masks for high contrast imaging. 
Astrophysical Journal, 590:593-603, June 10 2003a. 

R.J. Vanderbei, D.N. Spergel, and N.J. Kasdin. Circularly symmetric apodization via star- 
shaped masks. Astrophysical Journal, 599(1) :686-694, December 10 2003b. 

N. J. Kasdin, R. J. Vanderbei, M. G. Littman, and D. N. Spergel. Optimal one-dimensional 
apodizations and shaped pupils for planet finding coronagraphy. Applied Optics, 44 
(7):1117-1128, 2005. 

R. J. Vanderbei, N. J. Kasdin, and D. N. Spergel. Rectangular-mask coronagraphs for 
high-contrast imaging. Astrophysical Journal, 615, November 1 2004. 



R. A. Brown. Photometric and obscurational completeness. Astrophysical Journal, 624: 
1010-1024, 2005. 



R. A. Brown and C.J. Burrows. On the feasibility of detecting extrasolar planets by reflected 
starlight using the hubble space telescope. Icarus, 87:484-497, 1990. 

C. J. Burrows. Sharpness. In P. Y. Bely, editor, The Design and Construction of Large 
Optical Telescopes, pages 433-435. Springer- Verlag, New York, 2003. 

R.G. Lyon, J.M. Hollis, and J.E. Dorband. A maximum entropy method with a priori 
maximum likelihood constraints. Astrophysical Journal, 478:658-662, 1997. 



In this appendix we briefly summarize the derivation of Eq. 15. This is a slightly 
expanded treatment of the same problem in Burrows (2003). It is also demonstrates the 
overdetermined, MLE case in Lyon et al. (1997). As presented there, the problem is to 
estimate the irradiance, or photon count, of a signal corresponding to a point source in a 
known, uniform background. The image is given by the measurement Iij — B + AP^ + N^, 
where B is the assumed known background, A is the source irradiance, Py is the point spread 
function spread over the pixels and JV^- are noise photons in each pixel with zero 

mean and variance afy The noise is assumed independent and identically distributed across 
the pixels, so that E{NijN k i} = a^SikSji. What makes this problem more complicated than 
the standard linear least squares estimate is that the noise is Poisson distributed, so the 
variance at pixel is given by of - = B + APij, (where we are ignoring the readout noise). 



Since the background B is known, the weighted least squares estimate is found by 
minimizing the following cost with respect to A, 



A. Photometry via PSF Fitting 



A.l. Linear Least Squares Estimate 




(Al) 



where the Wij are some to be determined weights. 



This preprint was prepared with the AAS IATeX macros v5.2. 
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This results in the optimal estimate, 

, _ E„fe - B)W ijPij 
- w a p ij (A2) 

The optimal weights are found by forming the variance of the estimate and minimizing with 
respect to Wij (minimum variance estimation), 



r2 Y^mm 

A (En m 3 p*y 



(A3) 



The result is that oc 1/cr^-. That is, as is well known, the optimal weights are given by 
the noise variance. This implies that the cost function to be minimized is given by, 

ij j 

This equation has theoretical appeal since for large photon counts, the Poisson noise ap- 
proaches a Gaussian distribution and this is related to the maximum likelihood estimate. 
Note, however, that A appears nonlinearly in the cost and an exact maximum likelihood 
solution requires a more complex calculation. We address that further below. 

It is important to note that there is a bit of circular reasoning here. The expression for 
A in Eq. A2 was found by minimizing the cost with respect to A assuming the weights were 
independent of A. In the standard problem formulation, this is usually the case and the 
resulting estimate in Eq. A2 is optimal. In this case, however, this result was used to find 
that the optimal weights were proportional to the noise variance, which for Poisson noise is 
a function of A\ This puts into doubt the optimality of the resulting expression. However, 
the case of most interest here is the background limited case where B ^> APij. Here, the 
weights become constant and thus cancel out of the estimate. The resulting estimate is then 
given by, 

A= -' J ^ Tyl ' (AS) 



Zi jihj - B )P ij 

y.p 2 



and the variance of the error is given by, 

V P 2 cr 2 R 

2 _ £--nj V V _ D ( Kp\ 

a A ~ ( sp p2\2 ~ y p2 

where, in this background limited case, afj = B. 

A typical measure of the estimate quality is the signal-to-noise ratio given by dividing 
A by the standard deviation of its estimate and squaring, resulting in, for the background 
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limited case, 

Despite being non-optimal for the case where the background and planet are comparable, 
Eq. A5 has the appeal of being unbiased; it is therefore often used as the estimate for the 
general case as well. It is tempting to attempt a nonlinear least squares solution of the 
cost in Eq. A4, hoping that, though nonoptimal (the weights don't represent the minimum 
variance because the dependence on A was ignored), that it will still provide an improved 
SNR. Unfortunately, this is not the case. Not only is the problem much more difficult, with 
no simple closed form solution, the resulting estimate is biased and the variance is larger. 



A. 2. BLUE estimate 

In his appendix, Burrows (2003) used the BLUE (.Best Linear C/nbiased Estimate) ap- 
proach rather than least squares minimization. In fact, this is sensible as it avoids some 
of the circular reasoning problems above. We reproduce that here and show that the same 
signal-to-noise ratio results. 

The most general linear, unbiased estimate of the irradiance is given by, 



£ WijPi. 



(A8) 



It is straightforward to show that this is an unbiased estimate for any Wij as long as is 
zero mean. The variance of this estimate is found via, 



E{{A-Af} = 



E tj Eh hkhafWjjWki, 



V a 2 W 2 - 

°* = (STOP (A9) 

This is the same expression as in Eq. (D.l) of Burrows (2003). Following the same pro- 
cedure as above, the optimal weights are found by maximizing the signal-to-noise ratio or, 
equivalently, minimizing the variance. The resulting expression for the optimal weights is, 

(AlO) 

This shows that the weights are proportional to Pij/af-, as the term in parantheses is inde- 
pendent of Since any constant term cancels out of the estimate, the weights can be 
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set equal to Pij/crfj without loss of generality, 

w " = a^Tb < A11 > 

The circularity problem of Poisson processes shows up here when this expression for the 
weights is substituted into the BLUE estimate in Eq. A8, as the expression for A depends 
upon A itself. This is because the exact problem is inherently nonlinear; the linear estimate 
is only unbiased if the weights are forced to be independent. For an exact optimal estimate, 
one must turn to maximum likelihood, or other, approaches. However, again assuming a 
background limited signal, where APy <C B, the weights can be taken to be Wij = Pij. 
Substituting this into Eq. A8 results in the estimate equation, 

A = Ull £ p f )PiJ (M2) 

This is exactly the same as Eq. A5 for the least squares estimate and what we use for the 
PSF fitting detection in Eq. 15. 

As before, the expression for the weights, = P^, can be substituted into the variance 
expression (Eq. A9 above) to form a signal-to-noise ratio equation, 

It is interesting to look at the opposite case where the background is negligible. In 
that case, = 1 /A, a constant that cancels out of Eq. A8. This is equivalent to the least 
squares solution in the previous subsection weighted only by the PSF. This is a very common 
photometric approach to PSF fitting. The estimate simply becomes, 

A = ^ (A14) 

with variance, 

«a = v4r ( A15 ) 



This is an unbiased estimate with SNR of ^/ Aj^Pij- 

It is important to note that the estimator in Eq. A12 is still unbiased even if the 
background and signal are comparable; it is just no longer minimum variance. It is for that 
reason that it is used as the estimator for the detection problem. It interesting, however, 
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to compute the variance and resulting signal-to-noise ratio of the estimator in Eq. A12 but 
with the full noise variance. The expression for the estimate variance becomes: 

< = 45tfl + A (A16) 

resulting in an expression for the signal-to-noise ratio including the point source photon 
noise, 

(AV A'G^FZr 



Because of the Poisson nature of the signal and the resulting dependence of the noise 
variance on the signal itself, the linear, unbiased estimator is no longer optimal. Nevertheless, 
the linearity of the equation makes performing a detection and finding integration time 
estimates straightforward. It is, therefore, the most common approach to PSF fitting. It 
is, however, appealing to investigate the possibility of a maximum likelihood estimator that 
accounts for this dependence. We summarize such an approach next. 



A. 3. Poisson Maximum Likelihood Estimation 

Since the photons are known to arrive via a Poisson process, it is sensible to turn to 
maximum likelihood estimation under the Poisson distribution in an attempt to solve the 
optimality problem. 3 For simplicity, we start by examining the zero background case. The 
Poisson probability distribution at pixel (i,j) is given by, 

Vij =e- AP *^p- (A18) 

The total measurement probability across all the pixels is the product of the individual 
probabilities, as the measurements at each pixel are independent: 

p = e - £ ^ n(_^)^ (A19) 

Taking the logarithm gives the log-likelihood function, 



3 It is tempting to first search for a ML estimator assuming a Gaussian distribution of the noise, particularly 
since for large photon counts the Poisson distribution is close to Gaussian. However, Gaussian maximum 
likelihood in this case is difficult, highly nonlinear, and results in a biased estimate with larger error than 
the approximate linear cases. 
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where we have dropped constant terms independent of A. Taking the derivative and setting 
equal to zero gives the maximum likelihood estimate, 

A = St (A21) 

This is identical to the BLUE estimate with zero background in Eq. A14 and the PSF 
weighted linear least squares. For signals with little to no background, this is the most 
common optimal estimate of the intensity. 

For the Poisson maximum likelihood with background, we form the total probability 
density function, 

p = e -ZAP ij+ NB n(APv + BY" (A22) 

where iV is the number of pixels used in the estimate sums. We again form the log-likelihood 
function, 

L = - J2 Ap ij + Yl hi H Ap ij + B ) ( A23 ) 
and take the derivative and set equal to zero to find the equation that A must satisfy, 

y JijPij =y P (A24) 

For the general case of comparable background and signal, there is no closed form 
solution to this expression for A. One must resort to numerics. However, taking the expected 
value of both sides does show that the estimate is unbiased. At the cost of a loss of linearity, 
this is the best PSF fitting estimate for A, though it is difficult to find an expression for 
the signal-to-noise ratio that can be used to determine integration times. Monte-Carlo 
approaches are the most convenient. 

For the background limited case (B 3> AP^) it is possible to find an approximate 
expression for A by expanding Eq. A24, 



B 



(j2lijPiJ-Bj2 P ij) (A25) 



Again, the assumption is that the background is known. Because of the first order expansion, 
this estimate has a bias of 0(1/ B), small for the background dominated case. Because this 
ML expression does not lead to useful expressions of the integration time, we continue to 
use the (perhaps) more conservative linear approach. 
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Fig. 3. — Histograms of 50,000 Poisson simulations of the Bayesian test statistic for each 
case, pixels without a planet (Left) and pixels with a planet (Right). Dotted line indicates 
the threshold test value, 7. 
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Simulated SNR Test Statistic with and without Planet 
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Fig. 4. — Samples of the simulated PSF fitting test statistic with and without a planet for 
the selected integration time. 



